## Replication of
## Figure 3: Own Price Elasticities for Test Coffees and Competitor Coffees
## Figure B.1: Own Price Elasticities from AIDS Model

rm(list=ls())
library(foreign)
library(ggplot2)

# load data
dat <- read.dta("repFigure3boot.dta")

# logit model
m1 <- lm(delta~salesp+as.factor(storeprod)+time+time2,data=dat)
# elas
dat$elas <- coef(m1)["salesp"]*dat$salesshare*(1-dat$salesshare)*dat$salesp/dat$salesshare
aggregate(dat$elas,list(dat$prod),mean)

# block bootstrap logit model
dat$bstore <- NA
sims <- 5000
store <- matrix(NA,length(unique(dat$itemno)),sims)
n <- 26
set.seed(1234)
for(i in 1:sims){
print(i)
# sample stores
bo <- sample(unique(dat$store),n,replace=T)
# create boot sample
bdat <- dat[-c(1:nrow(dat)),]
for(p in 1:length(bo)){
btemp <-  dat[dat$store %in% bo[p],]
# boot store indicator
btemp$bstore <- p
bdat <- rbind(bdat,btemp)
}
# product store fixed effects
bdat$storeprod2 <- paste(bdat$itemno, bdat$bstore, sep = "-")
# fit model
m2 <- lm(delta~salesp+as.factor(storeprod2)+time+time2,data=bdat)
bdat$elas <-coef(m2)["salesp"]*bdat$salesshare*(1-bdat$salesshare)*bdat$salesp/bdat$salesshare
store[,i] <- aggregate(bdat$elas,list(bdat$prod),mean)[,"x"]
}
rownames(store) <- aggregate(bdat$elas,list(bdat$prod),mean)[,"Group.1"]

# block bootstrapped elas
bout <-rbind(apply(store,1,mean),apply(store,1,quantile,p=c(.05,.95)))
bout <- t(bout)
colnames(bout) <- c("el","lb","ub")
bout

## Figure 3: Own Price Elasticities for Test Coffees and Competitor Coffees 
dds1 <- data.frame(bout)
dds1$prod <- rownames(dds1)
dds1$plotn <- NA

dds1$plotn[dds1$prod=="COF FRENCH EXTRA DARK BULK"]  <- "FR Extra Dark"
dds1$plotn[dds1$prod=="COF COLOMBIAN SUPREMO BULK"]  <- "Colombian Supremo"
dds1$plotn[dds1$prod=="COF FRENCH ROAST OG BULK"]    <- "FR Regular "
dds1$plotn[dds1$prod=="COF BREAKFAST BLEND OG BULK"] <- "Breakfast Blend"
dds1$plotn[dds1$prod=="COF REGIONAL BLEND OG BULK"]  <- "Regional Blend"
dds1$plotn[dds1$prod=="COFFEE BLEND BULK"] <- "Coffee Blend "
dds1$plotn[dds1$prod=="COF MEXICAN OG BULK"]         <- "Mexican"

# elasticities from experiment
dds <- read.table("expelas.txt")
dds$prod  <- rownames(dds)
dds$plotn <- c("FR Regular","Coffee Blend")
dds <- dds[,-2]

# combine and plot
dds <- rbind(dds,dds1)
dds$plotn <- factor(dds$plotn,levels=c("FR Regular","Coffee Blend",
                                       "FR Regular ","Coffee Blend ",
                                       "Breakfast Blend","Colombian Supremo",
                                       "FR Extra Dark","Regional Blend",
                                       "Mexican")[9:1])
dds$gr <- c(1,1,rep(0,7))
dds$gr <- factor(dds$gr, levels=c(1,0),labels=c("Price Experiment","Historical Data"))

theme_bw1 <- function(base_size = 12, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) #% +replace%
  theme(
    axis.text         = element_text(size = base_size , colour = "black"),
    axis.ticks        = element_line(colour = "grey50")
  )
}

p = ggplot(dds,aes(x=plotn,y=el,fill=gr,shape=gr))
p = p + geom_point(size=4) + scale_y_continuous("Own Price Elasticity")
p = p +  scale_x_discrete("")  
p = p + geom_errorbar(aes(ymin=lb,ymax=ub,width=.4)) 
p = p + scale_shape_manual(name="Estimate from:",values=c(19,22))  
p = p + scale_fill_manual(name="",values=c("black","white"),guide="none")
p = p + coord_flip() + geom_hline(yintercept = 0,size=1,colour="black",linetype="dotted") 
p = p +  theme_bw1()
print(p)

dev.off()
pdf("fig3.pdf",width=8,height=7,pointsize=9)
print(p)
dev.off()


## replication code for AIDS Model
## Figure B.1: Own Price Elasticities from AIDS Model 
rm(list=ls())
library(foreign)
library(micEconAids)
library(ggplot2)

dat <- na.omit(read.dta("repFigureB1boot.dta"))

# functions to get elas
getela <- function(obj,set){
  if(set==1){
    out <- data.frame(rownames(obj$marshall)
                      ,diag(obj$marshall),diag(obj$marshallStEr)
                      ,diag(obj$hicks),diag(obj$hicksStEr)
    )
  } else {
    out <- data.frame(rownames(obj$marshall)
                      ,diag(obj$marshall),rep(NA,length(rownames(obj$marshall)))
                      ,diag(obj$hicks),rep(NA,length(rownames(obj$marshall)))
    )
    
  }
  return(out)
}

Shares <-  colnames(dat)[3:9]
Prices <-  colnames(dat)[10:16]
totExp <-  "totald"

# store dummies (-1), time, and time^2
Shifters <- c(colnames(dat)[17:41],"time","time2")

## fit AIDS
estResult <- aidsEst(Prices,Shares,totExp,method="IL",shifterNames=Shifters,
                     data = dat, maxiter = 1000 ,hom=T,sym=T)
## compute elasiticities
# get means for shifters
coefnew <- coef(estResult)
Zij <- as.matrix( apply(dat[,Shifters],2,mean))
coefnew[["alpha"]] <- as.vector(coef(estResult)$alpha+t(Zij)%*%t(coef(estResult)$delta))
names(coefnew[["alpha"]]) <- colnames(coef(estResult)$alpha+t(Zij)%*%t(coef(estResult)$delta))
coefnew[["delta"]] <- NULL
pMeans <- colMeans(dat[,Prices])
xtMean <- mean(dat[,totExp])
# elas
getela(aidsElas(coefnew, prices = pMeans, totExp = xtMean),set=0)[,c(1:2,4)]

## Blook Bootstrap 
sims <- 5000
storeh <- storem <- 
  matrix(NA,length(Shares),sims)
n <- 26
set.seed(1234)
dat$bstore <- NA
for(i in 1:sims){
  print(i)

# create boot sample
  bo <- sample(unique(dat$store),n,replace=T)
  bdat <- dat[-c(1:nrow(dat)),]
  for(p in 1:length(bo)){
    btemp <-  dat[dat$store %in% bo[p],]
    # boot store indicator
    btemp$bstore <- p
    bdat <- rbind(bdat,btemp)
    
  }
  bdat <- cbind(bdat,
                model.matrix(~as.factor(bdat$bstore))
  )
  
  Shifterd  <- colnames(bdat)[
    (which(colnames(bdat)=="(Intercept)")+1):ncol(bdat)
    ]
  Shifters <- c(Shifterd,"time","time2")
  
  ##  Fit AIDS
  estResult <- aidsEst(Prices,Shares,totExp,method="IL",shifterNames=Shifters,
                       data = bdat, maxiter = 1000 ,hom=T,sym=T)
  ## compute elasiticities
  # get means for shifters
  coefnew <- coef(estResult)
  Zij <- as.matrix( apply(bdat[,Shifters],2,mean))
  coefnew[["alpha"]] <- as.vector(coef(estResult)$alpha+t(Zij)%*%t(coef(estResult)$delta))
  names(coefnew[["alpha"]]) <- colnames(coef(estResult)$alpha+t(Zij)%*%t(coef(estResult)$delta))
  coefnew[["delta"]] <- NULL
  pMeans <- colMeans(bdat[,Prices])
  xtMean <- mean(bdat[,totExp])
  res3 <- getela(aidsElas(coefnew, prices = pMeans, totExp = xtMean),set=0)
  storem[,i] <- res3[,2]
  storeh[,i] <- res3[,4]
}

# plot results 
rownames(storeh) <- res3[,"rownames.obj.marshall."]

for(i in 1:length(rownames(storeh))){
  rownames(storeh)[i] <- 
    as.numeric(strsplit(rownames(storeh)[i],"esd")[[1]][2])
}

bout <- t(rbind(apply(storeh,1,mean),
                apply(storeh,1,quantile,p=c(.05,.95))))
colnames(bout) <- c("el","lb","ub")

dds1 <- data.frame(bout)
dds1$itemno <- as.numeric(rownames(dds1))
dds1$prod <- NA
dds1$plotn <- NA

dds1$prod[dds1$itemno==323]  <- "COF FRENCH EXTRA DARK BULK"
dds1$plotn[dds1$prod=="COF FRENCH EXTRA DARK BULK"]  <- "FR Extra Dark"

dds1$prod[dds1$itemno==325]  <- "COF COLOMBIAN SUPREMO BULK"
dds1$plotn[dds1$prod=="COF COLOMBIAN SUPREMO BULK"]  <- "Colombian Supremo"

dds1$prod[dds1$itemno==368]  <- "COF FRENCH ROAST OG BULK"
dds1$plotn[dds1$prod=="COF FRENCH ROAST OG BULK"]    <- "FR Regular "

dds1$prod[dds1$itemno==929]  <- "COF BREAKFAST BLEND OG BULK"
dds1$plotn[dds1$prod=="COF BREAKFAST BLEND OG BULK"] <- "Breakfast Blend"

dds1$prod[dds1$itemno==24177]  <- "COF REGIONAL BLEND OG BULK"
dds1$plotn[dds1$prod=="COF REGIONAL BLEND OG BULK"]  <- "Regional Blend"

dds1$prod[dds1$itemno==24189]  <- "COF WHOLE FOODS BLEND OG BULK"
dds1$plotn[dds1$prod=="COF WHOLE FOODS BLEND OG BULK"] <- "Coffee Blend "

dds1$prod[dds1$itemno==24191]  <- "COF MEXICAN OG BULK"
dds1$plotn[dds1$prod=="COF MEXICAN OG BULK"]         <- "Mexican"
dds1 <- dds1[,-4]


# elasticities from experiment
dds <- read.table("expelas.txt")
dds$prod  <- rownames(dds)
dds$plotn <- c("FR Regular","Coffee Blend")
dds <- dds[,-2]

# combine and plot
dds <- rbind(dds,dds1)
dds$plotn <- factor(dds$plotn,levels=c("FR Regular","Coffee Blend",
                                       "FR Regular ","Coffee Blend ",
                                       "Breakfast Blend","Colombian Supremo",
                                       "FR Extra Dark","Regional Blend",
                                       "Mexican")[9:1])
dds$gr <- c(1,1,rep(0,7))
dds$gr <- factor(dds$gr, levels=c(1,0),labels=c("Price Experiment","Historical Data"))

theme_bw1 <- function(base_size = 12, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) #% +replace%
  theme(
    axis.text         = element_text(size = base_size , colour = "black"),
    axis.ticks        = element_line(colour = "grey50")
  )
}

p = ggplot(dds,aes(x=plotn,y=el,fill=gr,shape=gr))
p = p + geom_point(size=4) + scale_y_continuous("Own Price Elasticity")
p = p +  scale_x_discrete("")  
p = p + geom_errorbar(aes(ymin=lb,ymax=ub,width=.4)) 
p = p + scale_shape_manual(name="Estimate from:",values=c(19,22))  
p = p + scale_fill_manual(name="",values=c("black","white"),guide="none")
p = p + coord_flip() + geom_hline(yintercept = 0,size=1,colour="black",linetype="dotted") 
p = p +  theme_bw1()
print(p)

dev.off()
pdf("figB1.pdf",width=8,height=7,pointsize=9)
print(p)
dev.off()






